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Abstract 

The internal energy of high-density hydrogen plasmas in the temperature 
range T = 10, 000 ... 50, OOOi^ is calculated by two different analytical ap- 
proximation schemes (method of effective ion-ion interaction potential - EIIP 
and Fade approach within the chemical picture - PACH) and compared with 
path integral Monte Carlo results. Reasonable agreement between the results 
obtained from the three independent calculations is found, the reasons for still 
existing differences is investigated. Interesting high density phenomena such 
as the formation of clusters and the onset of crystallization are discussed. 
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I. INTRODUCTION 



The thermodynamics of strongly correlated Fermi systems at high pressure are of grow- 
ing importance in many fields, including shock and laser plasmas, astrophysics, solids and 
nuclear matter, see Refs. for an overview. In particular, the thermodynamic properties 
of hot dense plasmas are essential for the description of plasmas generated by strong lasers 
Further, among the phenomena of current interest are the high-pressure compressibility 
of deuterium 0, metallization of hydrogen plasma phase transition etc., which occur in 
situations where both interaction and quantum effects are relevant. Among the early theo- 
retical papers on dense hydrogen we refer to Wigner/Huntington [§], Abrikosov Ashcroft 
10| and Brovman et al. |ll|] and, concerning the plasma phase transition, see Norman and 



Starostin [|T2| , Kremp et al. , Saumon and Chabrier ||T^ and Schlanges et al. |0 , as well 



as to some earlier investigations of one of us [p!6H T9[| . Among the early simulation approaches 
we refer to several Monte Carlo (MC) calculations, e.g. pO|- p^] . 

There has been significant progress in recent years in studying these systems analyt- 
ically and numerically, see e.g. [p],p|,^,p3|-p^ for an overview. However, there remains an 
urgent need to test analytical models by an independent numerical approach. Besides the 
molecular dynamics approach, e.g. p3| , |25| , the path integral Monte Carlo (PIMC) method 
is particularly well suited to describe thermodynamic properties in the region of high den- 
sity. This is because it starts from the fundamental plasma particles - electrons and ions, 
(physical picture) and treats all interactions, including bound state formation, rigorously 
and selfconsistently. We notice remarkable recent progress in applying these techniques to 
Fermi systems, for an overview see e.g. Refs. . 

Several methods have been developed to perform quantum MC. First we mention the 



restricted PIMC method [p0| -|33|; here special assumptions on the density operator p are 
introduced in order to reduce the sum over permutations to even (positive) contributions 
only. It can be shown, however, that this method does not reproduce the correct ideal Fermi 



gas limit p4[. An alternative are direct fermionic PIMC simulations which have occasionally 
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been attempted by various groups. Recently, three of us have proposed a new path integral 
representation for the N-particle density operator [p^ -p8| which allows for direct fermionic 
path integral Monte Carlo simulations of dense plasmas in a wide range of densities and 
temperatures. Using this concept, the pressure and energy of a degenerate strongly coupled 
hydrogen plasma have been computed PB|-P5[] as well as the pair distribution functions in 
the region of partial ionization and dissociation [^,^ . This scheme is rather efficient when 
the number of time slices (beads) in the path integral is less or equal 50 and was found to 
work well for temperatures ksT > O.lRy. 

One difficulty of PIMC simulations is that reliable error estimates are often not available, 
in particular for strongly coupled degenerate systems. Here, we will make a comparison 
with two independent analytical methods. The first is the method of an effective ion-ion 
interaction potential (EIIP) which has previously been developed for application to simple 
solid and liquid metals [pl] , p3| and which is here, for the first time adopted to dense hydrogen. 
The second is the method of Fade approximations in combination with Saha equations, i.e. 
the chemical picture (PACH) |^. The Pade formulas are constructed on the basis of the 
known analytical limits of low density and high density [Q, and they are exact up to 

quadratic terms in the density, interpolating between the virial expansions and the high- 
density asymptotics |T8|j4l|j4^ . 

We will show here that both methods, EIIP and PACH provide results for the internal 
energy which agree well with each other at high densities where the electrons are strongly 
degenerate and no bound states exist, approximately for n > lO^^cm^^. In this region, there 



is also good agreement with recent density functional results [^]. The agreement with the 
PIMC results is very good below lO^^cm"^. For intermediate densities, where the degree 
of ionization changes strongly, we observe deviations. Also, at high densities, the PIMC 
results, tend to lower energies than the analytical approaches. Finally, they reveal several 
interesting effects, such as the formation of clusters and the onset of ion crystallization. 
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II. PHYSICAL PARAMETERS AND BASIC EFFECTS 



Let us study a hydrogen plasma consisting of N^. electrons and Np protons {N^ = Np). 
The total proton (atom) density is n = Np/V. The average distance between the elec- 
trons is the Wigner-Seitz radius d = [3/47m]^/^, and other characteristic lengths are the 
Bohr radius ob = h/me'^, the Landau length / = e^/kT and the De Broglie wave length 
Ae = h/[2'KmckTY^'^ of the electrons. The degeneracy parameter is nAg. We define the 
dimensionless temperature r = kT/Ry which, in the considered below temperature interval, 
varies between 0.06 < r < 0.4. Furthermore, we introduce the Wigner-Seitz parameter 
Tg = d/as and the dimensionless classical coupling strength F = e^/{kTd). 

Hydrogen is anti-symmetric with respect to the charges (e_ = — e+) and symmetric with 
respect to the densities (n+ = n_ = n) and, due to the big mass difference, nip = 1836 me, 
ions and electrons behave quite differently. At the considered temperatures, the ions may be 
treated classically as long as n < 10^^ cm"'^. Further, for these temperatures and densities, 
the proton coupling parameter is in the range < F < 150, i.e. we expect strong coupling 
effects. We study in this work the internal energies of the fluid hydrogen system and start 
with providing some simple estimates for guidance. In the following we will give all energies 
in Rydberg units. 

First, at very low densities the electrons and the protons behave like an ideal Boltzmann 
gas. Therefore, the energy per proton (of free electrons and protons) is given by (in Rydberg 
units) 

e = E/N = 3r. (1) 

In other words the low-density limit is, in our temperature interval, a positive number in the 
region e ~ 0.2 — 1.2. With increasing density we expect a region where atoms and, possibly 
also a few molecules, are formed [|T^,0 . In the region of atoms a lower bound for the energy 
per proton is 



where the last term represents the binding energy IRy of H-atoms. If molecules are formed, 
the corresponding estimate per proton is still lower 

6= ^r-1.17. (3) 

Generally, the existence of a lower bound for the energy per proton was proven by Dyson 
and Lenard [4l| and Lieb and Thirring [45] 



E/N > -C, (4) 



where the best estimate known to us (which certainly is much too large), is C 23 [^]. We 
see that, with increasing density, the energy per proton tends to negative values and may 
reach a finite minimum. Further density increase will cause the energy to increase again as 
a result of quantum degeneracy effcts. 

In order to understand this increase let us look first at the limit of very high density 
(still in the region where the protons are classical). Then the first estimate of the energy is 

3 2.21 , . 

which is positive. The last term, representing the Fermi energy of the electrons, is strongly 
increasing with density (with power n^/^). In the next approximation according to Wigner's 
estimate we have to take into account the Hartree contribution to the electron energy and 
a corresponding estimate for the proton energy. The proton energy is estimated under the 
assumption that protons form a lattice. This way we find the estimate 



3 1.793\ /2.21 0.916 

r 



, . . , (6) 

2 rg J \ rj ^ ' 

The two corrections that were added to Eq. (|^) are both negative and scale like n^^^. In 
other words, these interaction terms might play a major role with decreasing density. At a 
critical density the energy per proton may become negative. This densitiy can be estimated 
from Eq. (|) by solving the quadratic equation 

= ^rr,2 - 2.709 r, + 2.21, (7) 



perturbatively, starting with the zero temperature hmit, and adding the first (hnear in r) 
correction, 

r° ~ 0.816 + 0.37r + .... (8) 

This result coincides, for r — 0, with Wigner's criterion for the existence of molecules: 
for d < as, molecules cannot exist since there is no room for forming bound state wave 
functions. According to Eq. (|^), for finite temperature, molecules exist only for still larger d 
as thermal fluctuations increase the wave function overlap. More generally, with increasing 
temperature, the energy becomes positive at lower density compared to the case T = 0. 

Summarizing the qualitative results obtained in this section we may state that we expect, 
in the given temperature range, the following general behavior of the internal energy per 
proton: at zero density the energy starts with the ideal gas expression which depends only 
on the temperature. With increasing density the energy per proton becomes negative due to 
correlation effects (bound states, electron correlations, proton correlations). A minimum is 
formed and at a density where the proton density is close to the inverse Bohr radius cubed 
the energy per proton turns to positive values and is more and more determined by the 
ideal electron energy increasing with n^/^, corrected by correlation contributions of order 
n^/'^ which are determined by the Hartree term and by proton-proton coupling effects. In 
the following we will show that this qualitative picture is supported by the results of our 
calculations. 

III. METHOD OF AN EFFECTIVE ION-ION INTERACTION POTENTIAL 

It is well known that in plasmas and plasma-like systems, in a broad parameter range, 
the interaction between the electron and ion subsystems is weak, whereas the interacrtions 
within the electron and ion subsystems can be strong. The corresponding small parameter 
is the ratio Uei/Ep of the characteristic value of electron-ion interaction Uei to the elec- 
tron Fermi energy Ep. Therefore, the mentioned approximation is valid for systems with 
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degenerate electrons, if Ep ^ > Ti, where Te and Tj are the electron and ion tempera- 
tures respectively (below we will consider the case Tg = Tj). Typical systems for which this 
approximation is fulfilled are simple solid and liquid metals and non-transitional metals in 
general, and this approximation serves as the basis for the computation of thermodynamic 
and electron kinetic properties, e.g. PB| , ^ |. 



For simple metals the Fermi energy is not very large compared to the characteristic 
electron-ion Coulomb interaction, taken at the average interparticle distance. However, due 
to the orthogonality of the wave functions for the conduction electrons and electrons bound 
in the ion shells, there is a partial compensation of the electron-ion Coulomb attraction at 
small distances which effectively weakens the electron-ion interaction. This fact is described 
in the theory of simple metals in the framework of the so called pseudopotential theory. The 
calculation of the pseudopotential is, in general, a complicated problem in particular due to 
its non-local structure |^6| , ^ . For practical applications it can be represented approximately 
as a local interaction with one or two fitting parameters for each metal. On basis of the 
pseudopotential theory all thermodynamic properties and electronic kinetic coefficients can 
be calculated with sufficiently high accuracy for a wide range of temperatures and pressures. 
Naturally, these calculations require reliable knowledge of the properties of the two quasi- 
independent subsystems: the degenerate electron liquid on the positive charge background 
and the classical ion subsystem with some effective strong inter-ion interaction. 

It is apparent that there is also a wide range of parameters for highly ionized strongly 
compressed hydrogen plasmas, where the electron-ion interaction is weak. For these param- 
eters the complicated problem of calculating the properties of a strongly coupled quantum 
electron-proton system can be essentially simplified. In so doing, the results obtained for 
high compression (when no bound electron states - hydrogen atoms and molecules - are 
existing), do not require any fitting, in contrast to the case of simple metals, because the 
inter-ion potential for hydrogen is pure Coulomb. Therefore, the data obtained with this 
analytical approximation, can be considered as an reliable basis for comparison with the re- 
sults of alternative approaches, including analytical and simulation methods for degenerate 



quantum systems of Fermi particles. The results of this pseudopotential approach are espe- 
cially important for conditions of extreme compression where the plasma is characterized by 
strong interaction within the electron and, especially, the ion subsystem. For these difficult 
situations experimental data are still missing whereas new acurate numerical methods for 
Fermi system are only emerging. 

Let us consider the Hamiltonian of an electron-proton plasma, for which the terms with 
infinite zero-components of the potentials are canceled, due to quasineutrality (for generality 
we retain the charge number Z of the ions): 

14^ 1 

k k,k',qjtO ^ k,q'jtO j=l 

Here is the energy of the electron with momentum hk and Uei{q) = is the Fourier- 

component of the electron-proton interaction potential. For the electron degrees of freedom 
in the Hamiltonian H the representation of second quantization is used where aj^ and Op are, 
respectively, the operators of creation and annihilation of an electron with momentum p . 
For the classical ions the coordinate representation is more convenient, thus in Eq. @ Ri 
denotes the coordinate of the i-th ion. To calculate the plasma energy, as in the theory of 



simple metals |11,23], two main approximations have to be used. The first is the adiabatic 



approximation for the ion motion, which is slow compared to the electron one. The second 

is the smallness of the ratio of the characteristic electron-proton Coulomb interaction to the 

Ze^ kT 

Fermi energy Ep. The respective parameter is V^i = - — = ZV — ~ n^^^^. Calculation of 

dEp Ep 

the electron energy in the external field of the immobile ions (protons) leads to the energy 
of the plasma given as function of the ion coordinates Rj. In general, the perturbation 
theory in terms of the parameter Tei gives rise not only to pair but, naturally, also to higher 
order ion-ion interactions, which are rather complicated. To second order of perturbation 
theory in the parameter Tei the energy per one electron of a plasma with a fixed proton 
configuration {Rj} is easily written, 
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+ ^EE^'^'(^)^"^'''"''^^- (10) 



2Ueiq = 0) 2VN, 

Here eg is the energy (per ion) of the correlated electron liquid on the homogeneous positive 
charge background. The functions Ile{q) and ee{q) are, respectively, the static polarization 
function and the static dielectric function of the correlated electron liquid. These functions 
are related to one another by the usual equality: 

ee(g) = l + ^ne(g). (11) 

The Fourier-component of the effective pair interaction potential between the ions, Vff , 
which appears in (|10D has the form: 

v-(,).i^-.i„)5^.m;!!. (12) 

r ^e{q) q ee[q) 

In the following, we will concentrate on hydrogen and set Z = 1 leading to the effective 
proton-proton interaction 

It is clear that, in contrast to liquid metals, where the presence of the pseudopotential 
leads to a more complicated structure of the effective potential, in a dense hydrogen plasma, 
the effective potential is determined only by electron screening. As it was shown in |[Tl[] for 
liquid metals, the additional pair interaction, arising from third and fourth order terms in 
the expansion of the electron energy in terms of the pseudopotential can play an important 
role in the effective interaction. For the effective potential of a hydrogen plasma a recent 
detailed analysis of these terms showed that these terms are essential only for rather 



rarified plasma conditions (r^ > 1.5), and they are practically negligible for higher densities, 
rs < 1.5, which we are considering in this paper. In fact, for > 1.6, the structure of 
the effective ion-ion potential in hydrogen changes drastically and can be considered as 
precursor of the appearence of molecular states. In this paper, we will use the simplest 
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version of the method of an effective ion-ion potential (EIIP) which includes the electron- 
proton interaction up to second order, so we are restricted to sufficiently high densities, 
corresponding to < 1.5. 

Further progress can be made by using for fig the random phase approximation (RPA), 
together with the long-wavelength and short-wavelength limits. 
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nRPA(^)^n^PA(o)M^ (14) 
o q 



where Hq^ = \/1mep is the Fermi momentum of the electrons. The analysis of this expression 
shows that the main contribution to the energy ([ToD comes from the small wave numbers. 
Therefore, with sufficient accuracy, it is possible to neglect the q-dependence of lie in Eq. (p!0D 
and, in particular, in the effective potential ([T^), replacing 11^^"^ (g) — > n^^'^(O). This means, 
we also neglect the well-known small oscillations of the effective potential for large distances, 
which are the result of a logarithmic singularity of the derivative [dX^^^^ j dq) \q=2qp- For the 
densities under consideration (which are much higher than usual metallic densities), these 
oscillations are not essential for the thermodynamic functions. At the same time, it is crucial 
to calculate the polarization function ne(0) fully selfconsistently: 

where eg is determined by ([ToD and, consequently, takes into account the electron-electron 
exchange and correlations. For the case of degenerate electrons we can use one of the 
analytical approximations for eg such as, for example, that of Nozieres and Pines or Wigner, 



see e.g. for an overview. Below we use Wigner's formula for the correlation energy, 
although for small the approximation of Nozieres and Pines is better (in fact, for the 
region < 1, where the deviations between these approximations for the correlation energy 
become essential, we can neglect correlations at all in comparison to kinetic and exchange 
terms). Because n^^'^(O) = K^pjiAa^e^^ it is clear that Eq. ( |T5| ) means renormalization of 
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-qrpa _^ -Q^ electron-electron interaction and, therefore, a renormalization of the 

momentum Ktf i^TF- 

, , /QTrV^^ 6 1 

Because for the considered approximation the effective proton-proton potential is described 
by the screened potential of Thomas- Fermi type, see Eqs. (|T^-(p!6D: 

2 

^pp{r) = —e~^, (17) 

we conclude that there is renormalization of the screening radius which is due to electronic 
correlations: 

rTF = - — = (18) 
Let us now rewrite Eg. (|10D for the considered approximation in the form: 

e = ee + ei (19) 

^. = lkT + — Y $™(i?, - R .) - £! f- + l-\ , (20) 

where k = d ■ ktf- After averaging over the proton positions with a Gibbs distribution 
(denoted by (...)), Eq. (p!9D can be represented as the sum of two terms: eg - the energy 
of a degenerate electron liquid on the positive homogeneous charge background and the 
energy of screened classical charged protons, interacting via the screened potential (0) and 
renormalized by the constant terms, obtained above: 

w + ?1 kT, (21) 



2 



with 



u 



(22) 
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Here u is the ionic interaction energy in kT— units. The energy (^) coincides with accuracy 
{kT/EpY with the usual thermodynamic energy determined from the free energy of the 
system because, in the considered parameter range, the electrons are degenerate (with the 
same accuracy). From expression ( pT]) follows that the energy of a classical one-component 
system of charged particles interacting via a screened (Debye or Yukawa) potential tends 
to infinity as SksTT /2k,'^ for k, (i.e. the screening radius diverges). The function u/T 
has been tabulated in |^|5T| (for the calculations of the phase diagram of a purely classical 
one-component Debye plasma), as function of the two parameters F and the dimensionless 
screening length k, based on accurate MD calculations for the Debye system. Below we use 
these numerical results to calculate the energy of a dense hydrogen plasma in the described 
above approximations. Within the Wigner approximation for the electron energy, 

'2.21 0.916 



+ ecorr Ry 



0.88 , , 

(23) 



7.8' 



we obtain, from Eq. ([Tq): 

22.1 



2 r 

_ 22.1 3.664 1.76r, 1.76r2 
'^^'''^ ~ VT ~ {rs + 7.8)2 - +7.«^3' ^24) 



where 7(rs — > 0) — * 1. Now, the total internal energy, Eq. (pT|) , can be expressed in terms 
of the tabulated function u/V as: 



Ry. (25) 



2.21 0.916 2 fu 3 

^1 ^ + '^°" + ^^f + 2f 

The numerical results computed from this approximation are included in Figs. |l|-0 below. 

Alternatively, we may use additional approximations for the computation of the internal 
energy of the plasma. This can be done by averaging Eq. ( PUj) over the ion Gibbs distribution 
with the same effective Hamiltonian (|TU]). Than we immediately find for the average energy 
per proton. 
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1 f (Pq 



where we introduced the ion- ion structure factor Sa^k) defined as 

^ ts I 1, A; = 

J [ 0, A; ^ 

Eq. ( P^D can be simphfied by replacing, approximately, the full structure factor by the 
OCP structure factor 5*^""^, computed with the effective ion-ion interaction. Then, the 
full energy can be written as the sum of three contributions: the first from the electron 
subsystem, the second from the classical ion OCP subsystem (both imbedded, respectively, 
into a positive and negative charge background) and a third term, ef*^^, which describes in 
perturbation-theoretical approximation for the polarization of the electron liquid by the ions. 
The resulting formulas coincide with the perturbation approximations derived by Hansen, 



DeWitt and others [21,221: 



^^™^=6e + 6P^^ + 5e, (28) 



oo 

o2 





As is clear from the above derivations, Eqs. (|28|), ( ^9]) are less accurate than the full EIIP 
model presented above. 

IV. FADE APPROXIMATIONS AND CHEMICAL PICTURE: PACH METHOD 



In this section we will explain in brief the method of Pade approximations in combination 
with the chemical picture, i.e. Saha equations P, p!8| , ^ , ^ (PACH). On the basis of the 
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PACH-approximation we will calculate the internal energy for the 3 isotherms T = 10, OOOK, 
30, OOOK, and 50, OOOK for those regions of the density where bound states (atoms and 
molecules) play a minor role. In other words we restrict our study to the density region 
where the plasma is strongly (but not necessarily fully) ionized. This method works only 
with analytical formulae which are, however, rather complicated; nevertheless the calculation 
of one energy data point takes no more than a few seconds on a PC. 

The Pade approximations were constructed in earlier work from the known analytical 
results for limiting cases of low density ||3|j40[| and high density The structure of the Pade 



approximations was devised in such a way that they are analytically exact up to quadratic 
terms in the density (up to the second virial coefficient) and interpolate between the virial 
expansions and the high-density asymptotic expressions |ll8|^^ . The formation of bound 



states was taken into account by using a chemical picture. 

We follow in large here this cited work, only the contribution of the OCP-ion-ion inter- 



action which is, in most cases, the largest one, was substantially improved following |^ 
With respect to the chemical picture we restricted ourselves to the region of strong ioniza- 
tion where the number of atoms is still relatively low and where no molecules are present. 
We will discuss here only the general structure of the Pade formulae. The internal energy 
density of the plasma is given by 

E = Eid + Eint (30) 

Here Eid is the internal energy of an ideal plasma consisting of Fermi electrons, classical 
protons and classical atoms and Eint is the interaction energy which is represented by 

Em = Np {te + e, + ta) (31) 

The splitting of the interaction contribution to the internal energy corresponds largely 
to the previous section. We have: 

• The electron-electron interaction: This term corresponds to the OCP energy of the 

we 



electron subsystem. Instead of the simple expressions used in earlier work ||18| , ^ , |39 
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used here a more refined formula for the energy . This formula is an interpolation 
between the Hartree limit with the Gellman-Brueckner correction (used already in the 
previous section), the Wigner limit and the Debye law including quantum corrections 



[r 



+ 50)r, + 2.3 + 2^6 di + rj 



(32) 



Here, a Wigner function has been introduced which is given by 

awix) = 2 6oa;log (l + [x^-^e-^'^/^^fco) + 2box/aw]~^) , (33) 

and the constants have the values d^ = 0.5; di = 0.6631; c?/^ = 0.125; an = 
0.91633; aw = 0.87553; bo = 0.06218; and hi = 0.0933. We mention that similar 
formulas are valid also for other thermodynamic functions by adjusting the constants 
The formula we have used here for the OOP contains all terms taken into account 
in the previous section but, in addition, also temperature dependent corrections. 

The ion contribution to the internal energy Cj. This term was calculated in the previous 
section. Here we will use a procedure which is based on the approximation (|27| , pS] ). 



This enables us to use results of the MC-calculations of Hansen, DeWitt and others 
2^,Q. According to Eqs. (p7| , ^81) the ion contribution is split into two terms 



e^ = + ef ^^ (34) 

where the first represents the O CP-contribution of the protons and the second the 
polarization of the proton OCP by the electron gas. For the region of high densities, i.e. 
large P and small we use the Livermore Monte Carlo data which were parametrized 
by DeWitt in the form 

^ocp ^ _,394gr + .8165P-^^ - .5012, (35) 
ef°^ = -Ts (.0543P + .1853P-^^ - .0659). (36) 
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We note that the polarization term describes the correction due to screening of the 
proton-proton interaction by the electron fluid. In order to obtain these expressions, 
semiclassical MC calculations were performed based on effective ion interactions which 
model the electrons as a responding background [^,^]. We do not need to go into 
the details of this method since the procedure corresponds to Eq. (^) derived in the 
last section. 

In the low density limit we used the Debye law with quantum corrections 

ef^^ = -.86603rrfor'-^[l - fiiF^'^], (37) 
gPOL _ _.7i744ri-5[i _ (38) 

Here, the temperature functions Bi and Ci describe rather complex quantum correc- 
tions which are, however, explicitly known and are easily programmed [0]. The Pade 
approximations which connect the high and the low density limits are constructed 
by standard methods |]l8|,p|,p2| and will not be given here in explicit form. For the 



O CP-energy of the ions we use the very accurate formulas proposed by Kahlbaum R 



The atomic contribution: In the region of densities and temperatures which is studied 
in this work this contribution gives only a small correction. We calculate the number 
of atoms on the basis of a nonideal Saha equation. The formation of molecules is 
not taken into account. We restrict the calculations to a region where the number 
density of atoms is so small, that the degree of ionization is larger than 75%. The 
contributions to the chemical potential which appear in the Saha equation are calcu- 
lated, in part, from scaling relations and, in part by numerical differentiation of the 



free energy given earlier |18,41]. For the partition function in the Saha equation we 



use the Brillouin-Planck-Larkin expression P,^. The nonideal Saha equation which 
determines the degree of ionization (the density of the atoms) is solved by 5-100 it- 
erations starting from the ideal Saha equation. Due to the high degree of ionization, 
the atomic interaction contributions can be approximated in the simplest way by the 
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second virial contribution and by treating the atoms as small hard spheres. 

The results of our Pade calculations for a broad density interval for three isotherms are 
included in Figs. |T]-^ 



V. SUMMARY OF THE PATH INTEGRAL MONTE CARLO SIMULATIONS 

The analytical approximations discussed in the previous sections work very well at high 
densities and if bound states are of minor importance. These conditions are not fulfilled for 
densities below than the Mott point corresponding to > 1. Here, recently developed path 
integral Monte Carlo simulations can be used. Starting from the basic plasma particles, 
electrons and ions, they "automatically" account for bound state formation and ionization 
and dissociation. Furthermore, in contrast to a chemical picture, no restrictions on the type 
of chemical species are made and the appearance of complex aggregates such as molecular 
ions or clusters of several atoms are fully included. On the other hand, the simulations are 
becoming increasingly difficult at high density where the electron degenercy is large. For 
this reason it is of high interest to compare results of the PIMC approach with alternative 
theories as they are expected to complement each other. This will be done in the next 
section. 

But first, we briefiy outline the idea of our direct PIMC scheme. All thermodynamic 
properties of a two-component plasma are defined by the partition function Z which, for the 
case of iVe electrons and Np protons, is given by 

Q(iVe,iVp,/3) 



Z(iV„iVp,V^,/?) 



with QiN„Np,P) = J2 j dqdrp{q,r,a;P), (39) 

" V 

where /? = l/ksT. The exact density matrix is, for a quantum system, in general, not known 



but can be constructed using a path integral representation 



(2) _pin) 



V " V 
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where p« = p {R^'-^\ R^'^ ; A(3) = |e-^^^|i?(^)), whereas A/5 = + 1) and AA^ = 

27rh'^Ap /rua, a = p,e. H is the Hamilton operator, H = K + Uc, containing kinetic and 
potential energy contributions, K and Uc, respectively, with (Jc = tj'j! + (j^ + U^^ being the 
sum of the Coulomb potentials between protons (p), electrons (e) and electrons and protons 
(ep). Further, = (g«,r«) = {Rf,Rf), for i = 1, . . . n + 1, = (g,r) = {Rf\Rf^), 
and = and a' = a. This means, the particles are represented by fermionic 

loops with the coordinates (beads) [R] = [R^^^; R^^^; . . . ; R^'^^; R^"^^^^], where q and r denote 
the electron and proton coordinates, respectively. The spin gives rise to the spin part 
of the density matrix S, whereas exchange effects are accounted for by the permutation 
operator P, which acts on the electron coordinates and spin projections, and the sum over 
the permutations with parity Kp. In the fermionic case (minus sign), the sum contains NJ/2 
positive and negative terms leading to the notorious sign problem. Due to the large mass 
difference of electrons and ions, the exchange of the latter is not included. 

To compute thermodynamic functions, the logarithm of the partition function has to be 
differentiated with respect to thermodynamic variables. In particular, the internal energy 
E follows from Q by 

f3E = -f3dlnQ/dl3, (41) 
This leads to the following result (for details, cf. [0), 



V s=0 

A^p ^ 2 " r A CI 2 Np N, 



p<t i'p*l p=i t=i ^IV 



1 ddetirj 



ab \s 



w,th CJ,= <|1*>. C^, = ^. (42) 

and = A/3a[/?'$"P(|xy,/5')]/<9/?'|/3'=A/3 contains the electron-proton Kelbg potential 
cf. Eq. (|45| ) below. Here, (. . . | . . .) denotes the scalar product, and qpt, Vpt and Xpt are 
differences of two coordinate vectors: qpt = qp — qt, fpt = Vp — rt, Xpt = Vp — qt, r^^ = Vpt + i/pt, 
= Xpt + Hp and = — y', with = AAg ^^^-^ ^a'^''. Here we introduced dimen- 
sionless distances between neighboring vertices on the loop, thus, explicitly, 

[r] = [r; yi^^; yP'*; . . .]. Further, the density matrix ps in Eq. (^21) is given by 

n Ne 

Ps{q, [r],P) = e-Z^^^^'M-^) J] H ^^^t |tf U (43) 

where U{q, [r], p) = U^iq) + {U%[r], AP) + U'P{q, [r], A/3)}/(n + 1) and 4)^^ = exp[-7r|ei'^|2]. 
We underline that the density matrix ( ^3D does not contain an explicit sum over the per- 
mutations and thus no sum of terms with alternating sign. Instead, the whole exchange 
problem is contained in a single exchange matrix given by 

W^ab \\s=\\e II.. (44) 

As a result of the spin summation, the matrix carries a subscript s denoting the number of 
electrons having the same spin projection. 

The potential appearing in Eq. (H2]) is an effective quantum pair interaction between 



two charged particles immersed into a weakly degenerate plasma. It has been derived by 
Kelbg and co-workers ||58|j5^ who showed that it contains quantum effects exactly in first 



order in the coupling parameter F, 

$"''(|r,,|,A/3) = -^ |l-e-^- + y^x,,[l-erf(x,,)]|, (45) 

■^abXab ■' 

where Xab = {'"'abl/ ^ab, and we underline that the Kelbg potential is finite at zero distance. 



The structure of Eq. (^) is obvious: we have separated the classical ideal gas part 
(first term). The ideal quantum part in excess of the classical one and the correlation 
contributions are contained in the integral term, where the second line results from the 
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ionic correlations (first term) and the e-e and e-i interaction at the first vertex (second and 
third terms respectively). Thus, Eq. (^) contains the important limit of an ideal quantum 
plasma in a natural way. The third and fourth lines are due to the further electronic 
vertices and the explicit temperature dependence [in Eq. (|^)] and volume dependence (in 
the corresponding equation of state result) of the exchange matrix, respectively. The main 
advantage of Eq. (^) is that the explicit sum over permutations has been converted into 
the spin determinant which can be computed very efficiently using standard linear algebra 
methods. Furthermore, each of the sums in curly brackets in Eq. (|42|) is bounded as the 
number of vertices increases, n —>■ oo. The error of the total expression is of the order of 



1/n. Thus, expression (^2]) and the analogous result for the equation of state are well suited 



for numerical evaluation using standard Monte Carlo techniques, e.g. [20,28 



In our Monte Carlo scheme we used three types of steps, where either electron or proton 

(k) 

coordinates, or qi or inidividual electronic beads Q were moved until convergence of the 
calculated values was reached. Our procedure has been extensively tested. In particular, 
we found from comparison with the known analytical expressions for pressure and energy 
of an ideal Fermi gas that the Fermi statistics is very well reproduced [^. Further, we 
performed extensive tests for few-electron systems in a harmonic trap where, again, the 
analytically known limiting behavior (e.g. energies) is well reproduced [|60|j6l[[ . For the 
present simulations of dense hydrogen, we varied both the particle number and the number 
of time slices (beads). As a result of these tests, we found that to obtain convergent results 
for the thermodynamic properties of hydrogen in the density-temperature region of interest 
here, particle numbers = Np = 50 and beads numbers in the range of n = 6 ... 20 are an 



acceptable compromise between accuracy and computational effort [36-38 
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VI. NUMERICAL RESULTS. COMPARISON OF THE ANALYTICAL AND 

SIMULATION DATA 



Let us now come to the numerical results. We have computed the internal energy of dense 
hydrogen using the two analytical (EIIP and PACH) approaches and the PIMC simulations. 
The data are shown in Figs. for three temperatures, 10, OOOK, 30, OOOK and 50, OOOK, 
respectively. 

Consider first the general behavior which is clearly seen for the lowest temperature, cf. 
Fig. [T|.b. The overall trend is an increase of the energy with density which is particularly 
rapid at high densities due to electron degeneracy effects; this is clearly seen from the ideal 
plasma curve (dash-dotted line). The nonideal plasma results show a prominent deviation 
from this trend which is in full agreement with the discussion given in Section ||: the 
formation of an energy minimum (where the energy may become negative) at intermediate 
densities. Our calculations for a nonideal hydrogen plasma asymptotically approach the 
ideal curve both, at low density (ideal classical plasma) and at high density (ideal mixture 
of classical protons and quantum electrons). At intermediate densities, between lO^^cm"^ 
and lO^^cm"^, the nonideal plasma energy is significantly lower than the ideal energy which 
is due to strong correlations and formation of bound states. In particular, we see clearly 
that indeed, for the considered temperatures, the total energy reaches negative values. 

Let us now compare the results from the different methods. First, we see that the energy 
minimum is reproduced by all methods, but there are quantitative differences regarding its 
depth and width. The general observation made for all temperatures, cf. also Figs. ^ and 
^, is that the simulations yield a deeper minimum and shift of the energy increase towards 
higher densities. Before further analyzing these differences, we concentrate on the results 
of the analytical approaches. For all temperatures (including higher ones), the PACH and 
EIIP approaches coincide in the limit of high densities. This is an important test since both 
contain the ideal Fermi gas result as a limiting case for high degeneracy. The interesting 
result is that this agreement holds up to densities as low as n = lO^^'cm"'^. For still lower 
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densities, the EIIP method yields lower energies which are closer to the PIMC results. At 
these densities, atom and molecule formation is becoming important, and both analytical 
methods (in their present form) are becoming unreliable. (For this reason, the Pade curve 
in Fig. 0.a is discontinued below lO^^cm"^, and in Fig. |^ the uncertain region is indicated 
by the dotted curve.) 

It interesting compare to another theoretical approach based on density functional (DFT) 
calculations. Recently, Xu and Hansen published data for T = 10, OOOK and < 



1.5 which are also included in Fig. |T|. These calculations which also neglect bound state 
formation practically coincide with the PACH results. The good agreement of the three 
completely independent approaches - EIIP, PACH and DFT - is a strong indication that 
they are able to yield reliable results for a fully ionized macroscopic hydrogen plasma at 
high densities, < 1.5. 

Let us now turn to the comparison with the PIMC simulations. As noted above, the 
overall agreement of all methods is satisfactory in view of the strength of correlation and 
quantum effects. Nevertheless, we observe deviations of our PIMC results from all other 
data, in particular around the energy minimum. Our data for T = 10, OOOK are also lower 
than restricted PIMC results of Militzer et al. cf. Fig. |T].b, whereas we found excellent 
quantitative agreement between the two independent quantum Monte Carlo methods above 



T = 50, OOOK, see the point for T = 62, 500K in Fig. ^.a, (see also Ref. [^). The reason for 
the low energies observed in our PIMC simulations at T = 10, OOOK are finite size effects: 
the homogeneous plasma state is unstable in the density region of the energy minimum. 
An analysis of the electron-proton configurations reveals that the plasma gains energy by 
forming small droplets which is a direct indication for a first order phase transition as 



discussed in the Introduction. These effects begin to appear in the weakly ionized plasma 
and are not contained in the present variants of the PACH and EIIP methods although they 
have been analyzed before ||65|. It is interesting to note that Xu and Hansen |Q observed 



strong fluctuations in their density functional calculations below = 1.5 which strongly 
resembled precursors of a phase transition. To clearify this interesting issue more in detail 
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requires extensive simulations which are presently under way. 

For completeness, we mention further effects which tend to lower the total energy and 
which are neglected in the analytical approaches: increased electron polarization and non- 
additive terms in the efficient proton-proton interaction which were analyzed by Kagan and 
co-workers ITT . 



The next interesting feature of the PIMC simulations is the shift of the energy growth to 
higher density values compared to the analytical models. This tendency becomes stronger 
with increasing temperature, as can be seen in Figs. There is no reason to doubt that the 
analytical methods (in accord with the density functional results) yield the correct energy 
asymptotics of a macroscopic electron-proton plasma at very high densities. (Account of 
proton degeneracy effects which are not included would only further increase the energies.) 
As noted above, our direct PIMC simulations become increasingly difficult with growing 
electron degeneracy, so we expect the results to become less accurate for densities exceeding 
lO^^cm"^. 

However, the most important effect results again from the finite-size character of our 
simulations. To better understand the high-density results, we analyze in Fig. ^ the electron- 
electron (e-e), proton-proton (p-p) and electron-proton (e-p) pair distribution functions. 
These functions exhibit features typical for strongly correlated systems. The most prominent 
effect is seen in the p-p function which exhibits a periodic structure at T = 50, OOOK which 
is even more pronounced at T = 10, OOOK. This proton ordering is typical for a strongly 
correlated ion fluid which is near the crystallization temperature ||66|. Our simulations for 
still higher densities show the formation of an ionic lattice immersed into a delocalized sea 
of electrons, i.e. an ionic Wigner crystal as it is known to exist in high density objects 
such as White or Brown dwarf stars. Thus, qualitatively, the simulations show the correct 
behavior at high densities. But due to the small size of the simulations (only 50 electrons and 
protons are presently feasible), the results are much closer to those for small ionic clusters 
which are known to exhibit quite peculiar behavior, including strong size dependence of 
the energy, negative specific heat etc. Therefore, in order to obtain more accurate data 
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for the internal energy of a macroscopic two- component plasma at ultrahigh compression, 
a significant increase of the simulation size is desirable which should become feasible in the 
near future. 



VII. DISCUSSION 

This work is devoted to the investigation of the thermodynamic properties of hot dense 
plasmas in the temperature region between 10, 0000 and 50, OOOK. We presented a new 
theoretical approach to high-density plasmas which is based on the theory of an effective 
ion-ion potential (EIIP). This method is shown to be quite efficient for fully ionized strongly 
correlated plasmas above the Mott density. 

Furthermore, a detailed comparison of several theoretical approaches on one hand and 
simulations on the other, has been performed over a wide density range. The first include 
the analytical models EEIP and the PACH on one hand and recent density functional data 



of Xu and Hansen on the other hand. The second group of data includes several new 
data points based on direct path integral Monte Carlo simulations (PIMC) of a correlated 
proton-electron system with degenerate electrons. In addition, we compared with restricted 



PIMC data of Militzer et al. 33 



From this comparison we conclude that the three theoretical approaches are in very good 
agreement with each other for a fully ionized hydrogen plasma in the high density region 
where < 1. On the other hand, the two simulations agree with each other for temperatures 
above 50, OOOK although no RPIMC data for high densities are yet available to us. This 
agreement over a broad range of parameters is certainly remarkable since the plasma is far 
outside the perturbative regime: it is strongly correlated and the electrons are degenerate. 
Moreover, all considered methods are essentially independent. 

Finally, the comparison of our PIMC simulation results with the analytical data reveals 
an overall good agreement, although deviations are observed above n = lO^^cm"^. The 
simulation energy reaches a far deeper minimum and the energy increase due to electron 
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degeneracy appears at higher densities. The discrepancy at lower densities (below n = 
lO^'^cm"'^) was attributed to not adequate treatment of bound state effects in the analytical 
methods, whereas the deviations at higher density are most likely due to finite size effects 
encountered by the simulations. These lead to droplet formation at low temperature and 
for densities between n = lO^^cm"^ and n = lO^^cm"^ which are an indication for the 
plasma phase transition |Q. At high density, the simulations reveal ordering of protons 



into a strongly correlated fluid and onset of the formation of a proton Wigner crystal. These 
interesting physical effects in high pressure hydrogen are of relevance for many astrophysical 
systems, but also for many laboratory experiments, including ultracold degenerate trapped 
ions and laser plasmas. 
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FIGURES 
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FIG. 1. Internal Energy of hydrogen for T = 10, OOOi^, a) normalized to the energy of a 
noninteracting electron-proton system and b) in units of 2N Rydberg. The curves show results of 
PACH-calculations ("Fade"), the EIIF model, our Monte Carlo simulations ("DPIMC"), density 
functional theory ("DFT") || and restricted FIMC data ("RFIMC") of Militzer et al. [||. 
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FIG. 2. Internal Energy of hydrogen for T = 30, OOOK. Same notation as in Fig. 1. Dotted 
line indicates region of low degree of ionization where the Fade and EEIP results are less reliable. 
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FIG. 4. Electron-electron (ee), proton-proton (pp) and electron-proton (ep) pair distribution 
functions of hydrogen from the PIMC simulations at n = lO^^cm"^ for a temperature of 10,000K 
(upper figure) and 50,000K (lower figure). Note the different vertical scales. 
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